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1.  INTRODUCTION 


In  radiation  transport  theory,  Monte  Carlo  techniques  have  been 
developed  that  reproduce  the  physical  model  in  as  much  detail  as  is 
necessary.  As  such,  these  approaches  provide  a powerful  tool  in  the 
solution  of  complex,  neutron  and  gamma-ray  transport  problems.  However, 
applications  outside  the  area  of  radiation  transport  have  been  dis- 
couraging. In  the  past,  brute-force  stochastic  approaches  have  been 
attempted  in  the  solutions  of  finite  difference  equations.  The  extrava- 
gant costs  of  development  and  computation  of  such  approaches  have 
maligned  Monte  Carlo  techniques  in  general.  More  recently1,  the  use 
of  modified  Green's  functions  has  produced  an  efficient  Monte  Carlo 
method  for  solving  complex  three-dimensional  time-dependent  parabolic 
and  elliptic  partial  differential  equations. 

This  paper  reports  a new  Monte  Carlo  approach  to  the  Dirichlet 
problem2 » 3»4 . This  approach,  which  simulates  a surface-density 
equation,  can  treat  problems  with  arbitrary,  geometric  complexity. 

It  requires  no  compromises  with  the  physics  of  the  problem  and  provides 
a new  insight  into  efficient  Monte  Carlo  solutions  of  many  difficult 
problems  of  mathematical  physics.  For  the  sake  of  clarity,  the 
terminology  of  heat  transfer  will  be  used,  but  this  surface-density 
approach  is  also  applicable  to  fluid  flow,  stress  analysis,  and  many 
other  problems  in  physics  and  engineering. 


1.  E.S.  Troubetzkoy  et  al.,  "Solution  of  Time  Dependent  Diffusion 
Equation  in  Complex  Geometry  by  the  Monte  Carlo  Method",  BRL 
Interim  Technical  Report,  MAGI  7053  (Rough  Draft)  1975. 

2.  T.J.  Hoffman  and  N.E.  Banks,  "Monte  Carlo  Solution  to  the 
Dirichlet  Problem  with  the  Double  Layer  Potential  Density", 

Trans.  Am.  Nucl . Soc.,  1£,  136,  (1974). 

3.  N.E.  Banks  and  T.J.  Hoffman.  "Continuous  Monte  Carlo  Solution 
to  the  Dirichlet  Heat  Transfer  Problem  Using  the  Double  Layer 
Potential",  20th  Army  Mathematics  Conference,  Boston,  MA  (1974). 

4.  T.J.  Hoffman  and  N.E.  Banks,  "Monte  Carlo  Solution  to  the 
Dirichlet  Problem",  Transactions  of  the  ANL-NEACRP  International 
Meeting  on  Monte  Carlo  Applications,  Argonne  National  Laboratory, 
July  1-3,  1974. 
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Section  2 describes  the  mathematical  bases  of  the  approach  and 
presents  a comprehensive  review  of  previous  Monte  Carlo  approaches. 
Section  3 develops  and  illustrates  the  concept  of  the  surface-density 
technique.  The  mathematical  development  involves  manipulation  of  the 
integral  equation  for  the  density  of  the  potential  of  a double  layer. 
Since  the  mathematical  justification  for  this  approach  is  intricate, 
a pictorial  illustration  will  be  presented  to  aid  computer  implemen- 
tation. Section  4 demonstrates  the  method  by  application  to  three 
heat-transfer  problems  of  varying  degrees  of  geometric  complexity.  The 
final  section  presents  the  conclusions  of  this  work  and  recommendations 
for  further  study. 


2.  BACKGROUND 

2.1  Problem  Statement. 

The  Dirichlet  problem  is  an  old  problem  of  mathematics  that 
frequently  arises  in  physical  applications.  The  problem  is  to  find 
a function,  T(r),  that  is  harmonic  in  a given  region  R,  i.e., 

V2T (r)  = 0.0  in  R (1) 

and  assumes  prescribed  values  on  the  boundary,  S,  enclosing  R,  i.e., 

T(?)  (rg)  on  S.  (2) 

From  a heat-transfer  point  of  view,  the  problem  is  to  find  the  steady- 
state  temperature,  T,  at  a point  r within  a homogeneous  region,  given 
the  surface  temperature,  <)>. 

There  are  many  ways  to  solve  this  problem  - both  deterministic 
and  stochastic.  Examples  of  deterministic  methods  are  finite  difference, 
finite  element,  Fourier  analysis,  and  Green's  functions.  These  methods 
will  not  be  discussed  here.  Instead  we  shall  describe  the  stochastic, 
or  Monte  Carlo,  approaches  to  the  solution  of  the  Dirichlet  problem. 

For  the  sake  of  clarity  these  approaches  will  be  illustrated  for  two- 
dimensional  regions,  but  all  are  applicable  to  three-dimensional 
problems . 

2 . 2 Discrete  Monte  Carlo  Approach. 

Perhaps  the  most  widely  used  Monte  Carlo  approach  to  the  Dirichlet 
problem  divides  the  region,  R,  with  a grid  as  shown  in  Figure  1 
(see  Curtiss  (5)  for  a literature  review).  The  temperature  at  each 


5.  J.H.  Curtiss,  "Sampling  Methods  Applied  to  Differential  and 

Difference  Equations",  Seminar  on  Scientific  Computation,  Inter- 
national Business  Machines  Corporation  (1949). 
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node  of  a uniform  grid  is  the  average  of  the  temperatures  at  the  four 
adjacent  nodes.  The  Monte  Carlo  random  walk  consists  of  starting  a 
"walk"  at  the  point  where  the  temperature  is  desired,  i.e.,  r in 
Figure  1.  With  equal  probability*  one  of  the  four  adjacent  nodes 
(indicated  by  "X's"  in  Figure  1)  is  selected,  and  the  walker  moves  to 
this  new  node.  The  temperature  at  the  selected  node  is  the  contri- 
bution to  the  temperature  at  r. 


Figure  1.  Discrete  Monte  Carlo  Approach  to  the  Dirichlet  Problem. 
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Since  the  temperature  at  the  new  node  is  not  known,  the  random-selection 
technique  is  iterated  until  the  walker  arrives  at  a point  where  the 
temperature  is  known,  i.e.,  a boundary  point  (see  Equation  (2)).  When 
the  walker  encounters  one  of  the  circled  nodes  which  approximate  the 


* The  mesh  size  need  not  be  constant  throughout  the  grid;  it  was  done 
here  for  the  purpose  of  illustration.  If  the  grid  is  not  uniform, 
the  probability  of  the  walker  moving  to  each  of  the  adjacent  nodes 
is  no  longer  1/4. 
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boundary,  this  boundary  temperature  is  saved  as  a contribution  to  the 
temperature  at  r.  Then  a new  walk  (history)  is  initiated  at  r,  and  the 
process  is  repeated.  The  average  of  these  saved  boundary  temperatures 
from  many  histories  is  a Monte  Carlo  estimate  of  the  desired  temperature 
at  r. 

With  this  method,  often  referred  to  as  the  "drunkard's  walk",  the 
walker  is  restricted  to  move  only  in  the  four,  specified  directions 
from  any  given  point.  Therefore,  only  the  finite-difference  approxi- 
mation to  the  problem  is  obtained,  so  that  the  accuracy  of  this  discrete 
approach  depends  on  the  fineness  of  the  grid.  Use  of  a finer  grid  not 
only  increases  the  number  of  possible  paths  by  which  a walker  can  reach 
the  surface,  but  also  allows  for  a more  accurate  representation  of  the 
boundary.  For  example,  if  a mesh  size  as  crude  as  that  shown  in 
Figure  1 were  used,  the  boundary  approximated  by  the  small  circles 
would  not  be  a very  accurate  representation  of  this  problem.  Of  course 
a finer  grid  structure  requires  a smaller  step  size  and  therefore  more 
steps  to  reach  the  boundary.  Consequently  the  computing  time  will  be 
longer. 

2 . 3 Floating  Random-Walk  Approach. 


* 


A Monte  Carlo  method  that  allows  large  random  steps  in  the 
solution  to  the  Dirichlet  problem  is  called  the  "spherical  process" 
or  "floating  random  walk".  This  method  was  proposed  in  1956  by 
Brown6  and  applied  to  heat  transfer  problems  in  1965  by  Haji-Sheikh 
and  Sparrow7*6.  Unlike  the  discrete  drunkard's  walk  previously 
described,  this  method  solves  the  Dirichlet  problem  as  expressed  in 
Equations  (1)  and  (2)  rather  than  the  finite-difference  approximation 
to  it.  This  floating  random-walk  method  is  illustrated  in  Figure  2. 

This  method  consists  in  constructing  a circle  centered  about  the 
point,  r,  at  which  the  temperature  is  desired.  This  circle  must  lie 
entirely  within  the  region,  R,  and,  for  efficiency,  should  be  tangent 
to  the  nearest  boundary  point  from  r as  shown  in  Figure  2.  The  first 
step  in  this  process  is  the  random,  uniform  selection  of  a point  on 
the  circumference  of  this  circle,  say  at  position  rj.  Another  circle. 


) 1 

r 


6.  G.M.  Brown,  "Monte  Carlo  Methods",  Modern  Mathematics  for  the 
Engineer,  McGraw-Hill  Book  Co.,  Inc.,  New  York,  NY  (1956). 

7.  A.  Haji-Sheikh  and  E.M.  Sparrow,  "The  Solution  of  Heat  Conduction 
Problems  by  Probability  Methods",  J.  Heat  Transfer,  ASME,  89 
(1967). 


8.  A.  Haji-Sheikh  and  E.M.  Sparrow,  "The  Floating  Random  Walk  and 
Its  Application  to  Monte  Carlo  Solutions  of  Heat  Equations",  J. 
SIAM,  Appl . Math.,  14,  (1966). 
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centered  at  r1  and  tangent  to  the  nearest  point  on  the  surface,  is  then 
constructed.  A point  is  randomly  selected  on  this  circumference,  e.g., 
r2,  and  the  walker  is  moved  to  this  point.  The  walk  continues  until 
it  encounters  the  surface,  and  the  surface  temperature  at  the  point  of 
encounter  is  saved.  A Monte  Carlo  estimate  of  the  temperature  at  the 
point  r is  the  average  of  many  such  random  walks.  Of  course,  the 
surface  is  never  really  encountered,  because  the  probability  of 
selecting  the  tangency  point  is  zero.  This  problem  is  circumvented 
by  defining  a thin  layer  inside  the  boundary,  indicated  by  the  dashed 
line  in  Figure  2.  When  the  walker  enters  this  layer,  he  is  assumed 
to  be  at  the  boundary. 


Figure  2.  Floating  Random-Walk  Approach  to  the  Dirichlet  Problem. 


The  floating  random  walk  allows  large  random  steps  for  points 
far  removed  from  the  boundary,  but  as  the  walker  approaches  the 
boundary,  smaller  steps  are  taken.  For  complex  geometries,  however, 
considerable  computer  time  is  required  to  locate  the  nearest  boundary 
point  from  the  current  walker  position  in  order  to  construct  the  circle. 
The  artificial  boundary  layer  required  for  convergence  also  introduces 
an  approximation,  or  error,  into  the  calculation.  If  it  is  too  large, 
the  wrong  problem  is  solved;  if  it  is  too  small,  the  computing  time 
is  extravagant.  The  error  approaches  zero  as  the  thickness  of  the 
layer  approaches  zero;  however,  the  computation  time  approaches  infinity. 
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Shapes  other  than  circles  have  been  proposed  in  an  attempt  to 
better  fit  the  boundary9.  A floating  random-walk  technique  using 
rectangles  (rectangular  parallelepipeds  for  three-dimensional  problems) 
is  under  development10.  This  approach  has  proved  very  useful  for 
nonlinear  problems  and  has  been  shown  to  be  very  efficient  in  geometries 
whose  physical  boundaries  are  rectangular. 


3.  THE  MONTE  CARLO  SURFACE-DENSITY  TECHNIQUE 


3 . 1 Mathematical  Development . 


Unlike  the  Monte  Carlo  approaches  previously  described,  the 
surface  density  approach  requires  no  artificial  boundaries  or  geometry 
searches  for  the  closest  boundary  point,  allows  large  random  steps, 
and  is  a continuous  method.  The  surface-density  approach  transforms 
the  differential  formulation  of  the  Dirichlet  problem,  Equations  (1) 
and  (2),  into  an  integral  formulation  using  the  density  of  the  potential 
of  a double  layer,  y(r  ) as  follows:11’12 


ft  = /d?s  v 


Cr-rs),ns 

4-n|r-r  |3 
' s 1 


where  y(rg)  is  the  solution  of  the  following  integral  equation. 


u(rs)  = 2<Krs)  + 


* / 


(r  -r  ) -n 
P sJ  r 


'dr  u(r  ) | - - 1 3 

p p r -r  3 
s v v 1 p s 1 


In  the  above  equations  drs  and  drp  are  surface  elements  about  the 
surface  position  vectors  rs  and  rp  respectively.  The  symbol  np 
denotes  the  unit,  inward-directed  vector  that  is  normal  to  S at 
rp.  T(r)  is  the  desired  temperature;  4>(rs)  is  the  known  surface 
temperature. 


M.E.  Muller,  "Some  Continuous  Monte  Carlo  Methods  for  the 
Dirichlet  Problem",  Annals  of  Math.  Statistics,  27_ , (1956). 

10.  M.H.  Kalos  and  H.  Steinberg,  private  communications  at  the  ANL- 
NEACRP  International  Meeting  on  Monte  Carlo  Applications, 
Argonne  National  Laboratory,  July  1-3,  1974. 

11.  S.G.  Mikhlin,  Linear  Integral  Equations,  Hindustan  Publishing 
Corporation,  Delhi,  India  (1960). 


12.  W.V.  Lovitt,  Linear  Integral  Equations,  Dover  Publications,  Inc., 
New  York,  NY  (1950). 
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The  notation  can  be  simplified  by  writing  Equations  (3)  and  (4) 


y(rs)  = Q(rs)  ♦ / drp  w(rp)  K(rp,rs) 


T(r)  = j drg  p(rs)R(r ,rs) , (6) 

where  the  definitions  of  Qfr  ),  K(r  ,r  ) and  R(r,r  ) follow  directly 
from  Equations  (3)  and  (4) . S P S 

The  reformulated  problem.  Equations  (5)  and  (6)  may  appear  more 
complex  than  the  original  problem.  Equations  (1)  and  (2) . The  moti- 
vation for  this  reformulation,  or  integral  transformation,  is  that 
similar  kinds  of  equations  are  solved  in  Monte  Carlo  radiation  transport. 
The  equation  that  is  normally  solved  in  Monte  Carlo  radiation  transport 
is  called  the  emergent-particle-density  equation13.  It  is  a Fredholm 
integral  equation  of  the  second  kind.  Equation  (5)  for  y(rs)  is  also 
a Fredholm  integral  equation  of  the  second  kind.  Stochastic  algorithms 
have  been  devised  for  solving  this  type  of  equation14*15. 

Consider  the  following  particle  transport  analogy: 

1.  X (register  for  the  quantity  of  interest)  is  set  to  0.0. 

2.  Q(rs)drs  is  the  probability  that  an  initial  event  will  occur 
in  drs  about  rs.  Q(r$)  is  assumed  to  be  a normalized, 
probability  distribution  function  (p.d.f.).  Therefore, 

on  the  boundary  of  the  region  an  initial  event  site,  r$0, 
is  selected  from  the  given  function  2<J>(rs).  The  quantity 
R(r,rso)  is  added  to  X . 

3.  K(rp,rs)d?s  is  the  probability  that  an  event  at  rp  will  be 
followed  by  an  event  in  drs  about  rs . Therefore,  another 
event  will  occur  with  probability  K(fso,?s)d?s  = a. 


13.  T.W.  Armstrong  and  P.N.  Stevens,  "A  V Importance  Function  for 
the  Monte  Carlo  Calculation  of  the  Deep  Penetration  of  Gamma 
Rays",  Pergamon  Press,  New  York,  NY  (1969). 

14.  H.  Kahn,  "Applications  of  Monte  Carlo",  Part  II,  AECU-3259,  Rand 
Corporation,  Santa  Monica,  CA  (1956)  . 


15.  R.R.  Coveyou,  V.R.  Cain  and  K.J.  Yost,  "Adjoint  and  Importance 
in  Monte  Carlo  Applications",  NSE,  27  (1967)  . 
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The  site,  r ,,  of  the  next  event  is  selected  from  the  p.d.f. 
si  r 

K(rsQ,rs).  The  quantity  R(r,rgl)  is  added  to  X. 

4.  Step  3 is  repeated  until  the  particle  fails  to  survive  an 
event  (probability  = 1-a). 

5.  When  the  particle  is  killed,  the  site  of  another  source  event 
is  selected  from  Q(r  ).  This  process  is  repeated  many  times, 
say  m. 

6.  A Monte  Carlo  estimate  of  T(r)  is  X/m. 


Two  major  drawbacks  soon  appear  when  the  above  procedure  is 
attempted  with  Equations  (5)  and  (6): 


1. 


Sampling  from  the  distributions  in  Equation  (5)  is 
extremely  difficult.  Except  for  very  simple  regions. 


selection  from  — K(rp,rg)  is  virtually  impossible. 


2.  The  Neumann  series  to  approximate  a solution  of  Equation  (5) 
(i.e.,  tallying  the  contribution  at  each  event  as  described 
in  the  particle  transport  analogy)  diverges12. 


The  first  of  these  difficulties  can  be  circumvented  with  a formulation 
adjoint  to  that  of  Equations  (5)  and  (6),  i.e., 


TW  = fdvs  Q(?s)  u*(?s),  (7) 

P*(rg)  = R(r,rs)  + J K(rg,rp)y*(rp)drp.  (8) 


With  Equation  (8)  initial  events  are  selected  from  R(r,r  ),  i.e., 
the  probability  of  the  first  event  occurring  in  the  surface  element 
d?s  is 


P(fs)d?s 


(r  - r ) -n 

s s , - 
dr 

4iv  | r-rg  | 3 


(9) 


Equation  (9)  is  a p.d.f.  for  convex  regions.  We  define,  for  brevity,  a 
convex  region  as  a region  comprised  of  a convex  set  of  points.  The 
solid  angle  subtended  by  the  surface  element  drg  when  viewed  from  r 
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r 


& 


is,  by  definition16, 


(r-r ) -n 

dn  = - - dr  . 

s 


(10) 


r-r 


Since,  for  convex  regions,  the  probability  of  an  initial  event  in  solid 
angle  dfi  about  the  direction  Q must  equal  the  probability  of  the  initial 


event  in  dr  about  r , 
s s 


p(fl)dfl  = p(rg)drs. 


(11) 


Substitution  of  Equations  (9)  and  (10)  into  Equation  (11)  gives 


1 


P(fi)  = 47  • 


(12) 


Consequently,  the  site  of  the  initial  event  can  be  determined  by 
selecting  an  isotropic  direction  from  r.  The  point  where  this 
trajectory  strikes  the  surface  is  r 


so 


Equation  (8)  makes  the  selection  of  subsequent  event  sites  as 
simple  as  the  selection  of  initial  sites.  If  an  event  occurs  at  r 


the  site,  r , , of  the  next  event,  can  be  selected  from 
si 


so 


lK(rsl»rso)ld'sl 


(r  - r , ) -n  . 
so  si'  si 


dr 


2tt  r - r . 

1 so  si1 


si 


(13) 


1 


The  right-hand  side  (RHS)  of  Equation  (13)  is  just  times  the  solid 


angle  subtended  by  drs^  as  viewed  from  rSQ.  Therefore,  subsequent 
sites  can  be  obtained  by  selecting  an  isotropic  direction  into  the 
region  from  rso.  The  point  where  a ray  in  this  direction  intersects 
the  surface  is  rsj.  Note  that  the  function  in  Equation  (13)  is  a 
p.d.f.  because  the  total  solid  angle  subtended  by  the  surface  from  a 
point  on  the  surface  is  2r  for  convex  regions. 


16.  W.  Kaplan,  Advanced  Calculus,  Addison-Wesley  Publishing  Co.,  Inc. 
Reading,  MA  (1959) . 


A statistical  weight  adjustment 


K(yy 


= -l  I is  also 


i^vy 


required  when  selection  of  subsequent  events  is  made  with 
Equation  (13). 
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The  above  treatment  of  the  selection  of  initial  and  subsequent 
event  sites  is  valid  only  for  regions  that  are  convex. 

However,  even  if  the  region  of  interest  is  not  convex.  Equation  (9) 
is  a properly  normalized  distribution  function.  Figure  3 illustrates 
this  point. 


Figure  3.  Selection  of  Initial  Event  Site  Using  Solid  Angle. 


The  probability  of  selecting  rso  on  the  surface  element  labeled  "1" 
is.  Pj  = ^ ; the  "probability"  for  surface  element  2 is  - ; and 

the  probability  for  surface  element  3 is  p . The  total  probability 


for  the  selection  of  r _ in  Afi  is 


Therefore,  the  initial  event  site  for  non-convex  regions  can  be 
selected  in  the  following  manner: 

1.  Select  an  isotropic  direction. 

2.  Save  all  boundary  positions  encountered  by  a ray  from  r in 
the  selected  direction.  (These  boundary  positions  are 
indexed  sequentially  as  encountered;  there  is  an  odd  number 
of  boundary  positions  for  all  closed  surfaces.) 

3.  Select,  with  equal  probability,  the  initial  event  site  from 
the  saved  boundary  positions. 

4.  If  the  index  of  the  selected  boundary  position  is  even, 
set  the  statistical  weight  of  the  "walker"  negative; 
otherwise,  set  it  positive. 

Selection  of  subsequent  event  sites  is  performed  in  a similar  manner. 

Consider  now  the  following  particle  transport  analogy  to  obtain 
a Monte  Carlo  solution  to  Equations  (7)  and  (8)  for  convex  regions: 

1 . X is  set  to  0.0. 

2.  A particle  isotropically  emitted  from  r strikes  the  surface 

at  r . 
so 

3.  Q(rgo)  added  to  X. 

4.  The  particle  is  reflected  isotropically  and  strikes  the 
surface  again  at  r 

5.  Q(rsl)  is  subtracted  from  X. 

6.  The  particle  is  reflected  isotropically  and  strikes  the 
surface  again  at  r^. 

7.  Q(rs2)  added  to  X. 


ad  infinitum. 

The  logic  of  the  above  random  walk  (other  than  the  "ad  infinitum") 
is  easily  adaptable  to  current  Monte  Carlo  radiation-transport  codes. 
Particles  are  simply  reflected  isotropically  from  the  inner  surface 
of  the  region.  The  calculation  depicts  particles  trapped  in  an 
internal  void  that  is  surrounded  by  a perfectly  elastic,  isotropic, 
albedo  surface. 


The  non-termination  of  the  random  walk  represents  the  divergence 
of  the  Neumann  series  approximation  to  Equations  (7)  and  (8) . This 
series  can  be  written  as 

N 

T(r)  = lim  'S*'  f dr  Q(r  ) p (r  ) , (15) 

N-k»  ' Vs  s n s 
n=0 

where 


u*  (r  ) = R(r,r  ) 
o b 

\ ^ = -/%  Ci(V|K(vy 1 • (16) 

The  limit  in  Equation  (15)  does  not  exist,  as  can  be  illustrated  in 
heat  transfer.  The  temperature,  $(r  ),  on  the  surface  of  a homo- 
geneous solid  is  given,  and  the  temperature,  T(r),  at  a point,  r , 
within  the  solid  is  desired.  Terms  in  the  Neumann  series.  Equation  (15), 
can  be  calculated  with  Monte  Carlo  methods  as  follows: 

1.  A large  number  of  particles,  m,  are  emitted  isotropically 

from  r . The  temperatures  at  the  points  these  particles 
strike  the  surface  [ <fr(r  ,) , $(r  2) , . . . ,$(r  ) ] are  added, 

multiplied  by  2.0,  and  divided  By  m.  Thissquantity,  say 

T is  an  estimate  of  the  first  term  in  Equation  (15). 

2.  These  m particles  are  re-emitted  isotropically  from  their 
impact  points.  Again  the  surface  temperatures  at  the  new 
impact  points  are  added,  multiplied  by  2.0,  and  divided  by 
m to  obtain  T1 . The  first  two  terms  in  the  series  of 
Equation  (15)  are  Tq  and  -T  . 

3.  This  process  is  repeated  N times.  The  series  can  now  be 
approximated  as 

N 

T(?  ) " £ (-1)\  * (17) 

n=0 

4.  Eventually  the  m particles  will  attain  an  equilibrium 
(eigenfunction)  distribution,  independent  of  r . Then 
the  expectation  value  of  Tn  will  equal  the  expectation 
value  of  Tn+j,  say  T . Therefore,  the  series  of  Equation 
(17),  and  hence  Equation  (15),  will  have  oscillations  of 
t Tg;  i.e.,  the  limit  does  not  exist. 
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Although  the  Neumann  series  does  not  converge  for  Equations  (7) 
and  (8),  it  can  be  shown  that  solutions  to  these  equations  do  exist 
and  are  unique11*12.  The  Neumann  series  is  not  the  most  general 
solution  to  the  Fredholm  equation,  i.e.,  there  are  Fredholm  equations 
that  have  solutions,  Equation  (8)  being  an  example,  and  yet  these 
solutions  cannot  be  expressed  as  a Neumann  series.  A series  can  be 
generated  for  the  Fredholm  integral  equation  as  follows: 


T(r)  = 


E/ 


dfs  Q(*s)  VV* 


where 


»*QCrs)  = R(r,is)  - hfdr p R(r ,fp)  | *(?,., ?p)  | 
u*l  (?s)  - - Jsl/dfp  R(r.rp)|K(rs,fp)|  -/drp  Jdi^i ,rp)  X 

|K(ft,rp) | |K(rs,rt) |]  (19] 

= '/  d?p  f°r  n > ^ 


This  series  can  be  shown  to  converge  to  the  solution  of  Equation  (18)  by  the 
alternating  series  test.  Implementation  of  the  above  solution  in  radiation 

codes  requires  changes  in  the  logic  used  for  the  Neumann  solution  to 
terminate  histories. 

The  proposed  method,  as  outlined  above,  appears  at  first  glance 
to  be  useful  for  the  calculation  of  T(r)  at  only  one  point  at  a time. 
However,  T(r)  can  be  calculated  at  many  positions  simultaneously,  with 
only  a small  increase  in  computational  time,  by  using  statistical 
weights.  For  example,  if  T(r)  is  desired  at  rj  and  r2 , the  initial 
event  can  be  selected  by  starting  a particle  isotropically  from  rj . 
Contributions  for  r^  are  then  multiplied  (weighted)  by 

__  I (?2  - ?,)•",!  I (*i  - *SISI  _ (20) 

R(flffs)  I (r:  - rg) -ns]  ( |r2  - rj3] 


where  rg  is  the  site  of  the  initial  event.  This  weighting  is  justified 
because  once  a particle  strikes  the  surface  at  a given  point,  the 
random  walk  is  independent  of  its  origin.  This  weighting  is  a highly 
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useful  property  of  the  surface  density  technique.  With  Equation  (20) 
the  temperature  at  any  number  of  points  (temperature  distributions)  can 
be  calculated  from  one  set  of  Monte  Carlo  histories.  The  temperature 
differences  (or  gradients)  calculated  in  this  correlated  fashion  are 
more  meaningful  than  those  from  independent  histories. 

3.2  Pictorial  Demonstration. 

The  surface  density  method  is  illustrated  in  Figure  4. 


Figure  4.  Surface  Density  Approach  to  the  Dirichlet  Problem. 


A walker  starts  at  the  point  where  the  temperature  is  desired,  r, 
in  Figure  4.  It  leaves  this  point  isotropically  and  moves  directly 
to  the  surface  with  no  intermediate  steps.  The  temperature  at  the 
surface  point  encountered  by  the  walker  is  saved.  The  walker  then 
moves  from  this  surface  point,  isotropically-inward,  directly  to 
another  surface  point.  The  walk  is  continued  until  four  or  five 
surface  temperatures  are  collected.  These  temperatures  are  combined 
in  a series  (Equation  (18)),  and  this  combined  temperature  is  a 
Monte  Carlo  estimate  of  the  temperature  at  r 1 . Since  the  series  for  the 
Fredholm  equation  is  truncated,  i.e.,  only  four  or  five  surface  tempera- 
tures are  used,  a numerical  approximation  results.  The  truncation  error 
can  be  estimated  from  the  last  calculated  term  of  the  series.  The 
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temperature  at  r2  in  Figure  4 can  be  obtained  by  weighting  each 
history  from  i-j  using  Equation  (20) . 


4.  APPLICATIONS 


4.1  Two-Dimensional  Convex  Regions. 


The  Monte  Carlo  radiation  transport  code  MORSE17  was  altered  to 
provide  isotropic  reflection  at  albedo  surfaces,  to  calculate  the 
terms  of  the  series  (See  Equation  (18)),  and  to  weight  the  detector 
positions  (see  Equation  (20)).  Since  the  first  problem  is  two- 
dimensional,  it  requires  a slight  modification  to  Equations  (8)  and 
(20).  The  modified  two-dimensional  forms  are  listed  below  for 
reference: 

, ^sH  [ (V'p^s 
v (r  ) = + / L v (r  )dr  , (8a) 

2tt | r-r  | 2 J »|r  -r, |2  P P 


R(yrs)  = [ (r2-rs)-ns][  ly?sl2l 

R(?r?s)  I (r1-?s)^sH  \i2-is\2] 


(20a) 


We  require  the  temperature  at  various  points  within  a square 
situated  in  the  region  0 < x < 1,  0 < y < 1.  The  surface  temperature 
is  unity  on  the  side  x=0,  0 < y < 1 and  zero  on  the  other  boundaries, 
Separation  of  variables  provides  the  analytic  solution, 

cc  (21) 

4 V ^ _i_L  r ..M  — j r n-,  ..l 


T(x,y)  = - 
ir 


(2n+l)sinh(  (2n+l)ir] 


Temperatures  were  calculated  at  nine  (detector)  positions  with  the 
initial  (source)  position  determined  by  random  selection.  The  Monte 
Carlo  results,  along  with  the  corresponding  analytic  results,  are 
shown  in  Table  I.  This  run  of  10,000  histories  required  17  s of  CPU 
time  on  the  CDC  7600.  An  identical  calculation  for  only  one  detector 
required  12  s. 


E.A.  Straker,  N.R.  Bym  and  W.H.  Scott,  Jr.,  "The  MORSE  Code 
With  Combinatorial  Geometry",  SAI-72-511 -LJ , Science  Applications, 
Inc.,  (1972). 
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TABLE  I.  RESULTS  OF  TWO  DIMENSIONAL  CALCULATION 


(a) 


detector 

POSITION  (x,y) 

ANALYTIC 

MONTE  CARLO 

1 

(.25,  .25) 

0.43 

0.43  ± 0.03 

2 

(.25,  .50) 

0.54 

0.54  ± 0.03 

3 

(.25,  .75) 

0.43 

0.45  ± 0.03 

4 

(.50,  .25) 

0.18 

0.20  ± 0.02 

5 

(.50,  .50) 

0.25 

0.26  ± 0.02 

6 

(.50,  .75) 

0.18 

0.19  ± 0.03 

7 

(.75,  .25) 

0.07 

0.08  ± 0.02 

8 

(.75,  .50) 

0.10 

0.11  ± 0.02 

9 

(.75,  .75) 

0.07 

0.08  ± 0.02 

(a)  These 

Monte  Carlo  results  were 

obtained  in  one 

run  of  10,000 

histories  which  required  17  s of  CPU  time  on  the  CDC  7600; 
temperatures  for  the  nine  detectors  were  obtained  simultaneously 
by  the  weighting  technique  (Equation  20a) . 


4 . 2 Three-Dimensional  Convex  Regions. 

The  second  problem  treats  a finite  cylinder  with  height  and 
radius  both  equal  to  6 metres.  The  outer  surface  of  the  cylinder  is 
at  400°C  with  the  top  and  bottom  surfaces  at  100  and  300°C,  respec- 
tively. The  analytic  solution  is  given  by  the  following  expression18: 


i 

I 


T(r,z)  = 100  ^ 

n=l 


VrV 


(6an)sinh(6an) 


sinh  (6-z)an 


sinh  a z -i 


16  I [ (2n-l)irr/6]  (2n-l)irz 

— sin 

tt  (2n-l)I  [ 2n-l)ir]  6 


(22) 


where  a are  the  roots  of  the  Bessel  functions, 
n 


18.  H.S.  Cars law  and  J.C.  Jaeger,  Conduction  of  Heat  in  Solids, 
Oxford  at  the  Clarendon  Press,  2nd  ed.  New  York,  NY  (1959). 


> • 
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Temperatures  were  calculated  at  ten  positions  with  the  source  weighting 
of  Equation  (20) . The  analytical  and  Monte  Carlo  results  are  compared 
in  Figure  5. 

4.3  Three-Dimensional  Non-Convex  Regions. 

The  third  problem  treats  a region  with  both  concave  and  convex 
surfaces,  viz.,  a finite  hollow  cylinder  (Figure  6).  MORSE  was  modified  to 
generate  secondary  particles  at  the  surface  of  the  region.  These  secondary 
particles  were  given  positive  or  negative  statistical  weights  in 
accordance  with  the  algorithm  that  follows  Equation  (14).  The  Monte 
Carlo  results  from  MORSE  and  the  analytical  results  from  Carsiaw  and 
Jaeger  (is)  in  Table  II. 

TABLE  II.  RESULTS  FOR  FINITE  HOLLOW  CYLINDER  PROBLEM. 


R 

Z 

M.C.a 

Analytic 

5.0 

5.0 

134°C  ± 

11°C 

134°C 

5.0 

4.0 

141°C  ± 

7°C 

142°C 

5.5 

3.0 

123°C  ± 

5°C 

121°C 

5.0 

3.0 

141°C  ± 

4°C 

144°C 

4.5 

3.0 

164°C  ± 

6°C 

174°C 

aThese  Monte  Carlo  results  were  obtained  with  one  run  of  10,000 
histories  which  required  2.04  min  of  CPU  time  on  the  CDC  7600. 


Three  independent  calculations  were  made  for  this  problem.  The 
first  case,  for  a single  temperature,  required  1.80  min  on  the  CDC 
7600  ; the  second,  for  two  temperatures,  1 .87  min.  The  third  calculation, 
for  the  temperatures  at  the  five  positions  shown  in  Table  II,  required 
2.04  min  of  computer  time.  These  times  demonstrate  the  advantages  of 
the  weighting  technique.  Equation  (20). 
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Results  for  Solid  Cylinder  Calculation 


Geometry  and  Boundary  Conditions  for  Finite  Hollow  Cylinder 


5 . SUMMARY 

The  surface-density  technique  is  a continuous  Monte  Carlo  method 
for  the  solution  of  the  Dirichlet  heat  transfer  problem.  The  only 
other  continuous  Monte  Carlo  method  that  allows  large  random  steps  in 
the  solution  to  this  problem  is  the  floating  random-walk  method. 

The  major  advantages  of  the  surface-density  walk  over  the  floating 
random-walk  are: 

1.  The  floating  random  walk  required  very  small  random  steps 
near  the  boundary.  The  surface-density  walk  proceeds 
directly  to  the  surface. 

2.  Determination  of  the  closest  boundary  position  (required 

in  the  floating  walk  at  each  event  site)  is  extremely  diffi- 
cult. The  surface  density  walk  requires  only  the  surface 
intercept  of  a walker  moving  in  a given  direction  from  a 
given  point.  Geometry  searches  of  the  latter  type  are  less 
time  consuming  even  for  simple  geometries;  for  complex 
geometries  this  advantage  becomes  more  pronounced. 

3.  The  surface-density  method  does  away  with  the  artificial 
boundary  layer  and  its  attendant  errors  required  for 
convergence  of  the  floating  walk.  The  surface-density 
walk  does  introduce  a truncation  error  in  evaluating  the 
series,  but  this  error  can  be  estimated  in  the  calculation. 

4.  Current  Monte  Carlo  radiation  transport  codes  are  readily 
adaptable  to  the  surface-density  method  but  not  to  floating 
walk.  The  geometry  search  logic  described  above  for  the 
surface-density  approach  is  standard  in  all  Monte  Carlo 
radiation  transport  codes. 

5.  The  surface-density  walk  provides  temperature  distributions 
from  one  set  of  Monte  Carlo  histories  through  statistical 
weighting.  The  floating  random  walk  provides  the  temperature 
only  at  one  point  with  a given  set  of  histories. 


The  surface-density  method  does  have  its  limitations.  For  heat 
transfer  these  include: 

1.  The  surface  density  approach  is  currently  limited  to  a 
single  medium. 

2.  The  type  of  boundary  condition  determines  the  integral 
formulation,  i.e.,  Equations  (7)  and  (8).  Therefore 
different  types  of  boundary  conditions,  e.g.,  convective 
boundary  conditions,  require  the  simulation  of  different 
integral  equations. 
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3. 


The  surface-density  Monte  Carlo  approach  as  formulated  in 
this  paper  is  applicable  only  to  time- independent  problems. 


These  limitations  severely  restrict  the  applicability  of  the 
surface-density  Monte  Carlo  technique.  However,  these  limitations 
can  perhaps  be  overcome  with  further  development.  Since  the  surface- 
density  Monte  Carlo  approach  is  readily  incorporated  in  Monte  Carlo 
radiation  transport  codes,  we  feel  that  such  development  is  worth 
pursuing.  This  compatibility  is  highly  useful  for  the  many  problems 
that  require  calculations  of  both  heat  transfer  and  radiation 
transport.  If  this  compatibility  can  be  maintained,  then  criticality, 
shielding,  and  heat  transfer  can  be  calculated  with  one  geometric 
description  and  one  computer  code. 

This  method  should  be  applied  to  heat  transfer  and  boundary  layer 
analyses  in  problems  of  interior  ballistics.  The  aspect  of  a boundary 
layer  being  physically  a single  medium  transport  phenomenon  makes  this 
approach  extremely  attractive.  Also,  making  use  of  this  technique  would 
allow  the  boundary  layer  analysis  and  the  heat  transport  processes  to  be 
simulated  by  a single  computer  code  without  recourse  to  gross  simplifying 
assumptions  which  are  usually  made  in  the  solution  to  boundary  layer 
problems. 
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